%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% This program replicates Figure 7 in the paper "Optimal 
%%% Portfolio Choice with Fat Tails and Parameter Uncertainty",
%%% Journal of Financial and Quantitative Analysis (2024),
%%% Raymond Kan & Nathan Lassance 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Note: due to the use of a finite number of simulations 
%%% to evaluate (k1,k2,k3), the final results are random 
%%% and will not correspond exactly to those in the paper.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear
close all

%%% Load data
load Dataset25SBTM.txt;
R = Dataset25SBTM;

%%% Define variables
nobs = 10^1;%We use nobs=10^3 to create Figure 7 in the paper
N = length(R(1,:));
Ttotal = length(R(:,1));

%%% Exponential integral
Expfun = @(n,x) integral(@(t)exp(-x.*t)./(t.^n),1,inf);

%%% T = 60
T = 60;
a1 = (T-N-2)/N;
a2 = a1*(T-N-1)*(T-N-4)/(T-2);
a3 = a2/T;
M = eye(T)-(ones(T,1)*ones(1,T))./T;
rho = N/T;
NumberWindow = Ttotal-T;
EtaHat = zeros(NumberWindow,1);
VarphiHat = zeros(NumberWindow,1);
k1Hat = zeros(NumberWindow,1);
k2Hat = zeros(NumberWindow,1);
k3Hat = zeros(NumberWindow,1);
EtaTilde = zeros(NumberWindow,1);
VarphiTilde = zeros(NumberWindow,1);
k1Tilde = zeros(NumberWindow,1);
k2Tilde = zeros(NumberWindow,1);
k3Tilde = zeros(NumberWindow,1);
for j = 1:NumberWindow
    if mod(NumberWindow-j,100) == 0
        WindowsLeftT60 = NumberWindow-j
    end    
    X = R(j:j+T-1,:);
    %(\eta,\varphi) Student t
    [~,~,nuhat] = mlstudent(X);
    nuhat = min(nuhat,1000);
    y = @(x) (nuhat-2).*rho.*x./(2*(1-rho));
    fun = @(x) y(x).*exp(y(x)).*Expfun(nuhat/2,y(x))-rho;
    EtaHat(j,1) = fzero(fun,0.5*(1+nuhat/(nuhat-2)));
    VarphiHat(j,1) = 2*EtaHat(j,1)^2*(1-rho)/(nuhat-EtaHat(j,1)*(nuhat-2));
    %(\eta,\varphi) Elliptical
    muhat = mean(X);
    sigmahat = cov(X,1);
    tauhat = sum((X-muhat).^2,2);
    tauhat = tauhat./mean(tauhat);
    fun = @(x) sum(1./(T-N+x.*tauhat.*N))-1;
    EtaTilde(j,1) = fzero(fun,[1,10]);
    VarphiTilde(j,1) = (1-rho)./(1/EtaTilde(j,1)^2-sum((N.*tauhat.^2)./...
                       (T-N+N.*EtaTilde(j,1).*tauhat).^2));
    %(k1,k2,k3) Student t
    Expk1vec = zeros(nobs,1);
    Expk2vec = zeros(nobs,1);
    Expk3vec = zeros(nobs,1);
    for n = 1:nobs
         lambdahat = diag(sqrt((nuhat-2)./chi2rnd(nuhat,[T,1])));
         B = lambdahat*M*lambdahat;
         oneslambda = lambdahat*ones(T,1);
         Y = randn(T,N);
         Mat1 = inv(Y'*B*Y);
         Mat2 = Mat1^2;
         Expk1vec(n,1) = trace(Mat1);
         Expk2vec(n,1) = trace(Mat2);
         Expk3vec(n,1) = oneslambda'*Y*Mat2*Y'*oneslambda;
    end
    k1Hat(j,1) = a1*mean(Expk1vec);
    k2Hat(j,1) = a2*mean(Expk2vec);
    k3Hat(j,1) = a3*mean(Expk3vec);
    %(k1,k2,k3) Elliptical
    lambdahat = diag(sqrt(tauhat));
    B = lambdahat*M*lambdahat;
    oneslambda = lambdahat*ones(T,1);
    Expk1vec = zeros(nobs,1);
    Expk2vec = zeros(nobs,1);
    Expk3vec = zeros(nobs,1);
    for n = 1:nobs
         Y = randn(T,N);
         Mat1 = inv(Y'*B*Y);
         Mat2 = Mat1^2;
         Expk1vec(n,1) = trace(Mat1);
         Expk2vec(n,1) = trace(Mat2);
         Expk3vec(n,1) = oneslambda'*Y*Mat2*Y'*oneslambda;
    end
    k1Tilde(j,1) = a1*mean(Expk1vec);
    k2Tilde(j,1) = a2*mean(Expk2vec);
    k3Tilde(j,1) = a3*mean(Expk3vec);
end
fig = figure;
set(groot,'defaultAxesTickLabelInterpreter','latex');  
subplot(3,3,1)
boxplot([EtaHat,EtaTilde,VarphiHat,VarphiTilde],'Symbol','');
set(gca,'TickLabelInterpreter','latex');
grid on
title('$T=60$','interpreter','latex');
ylim([0.95 12]);
xticklabels({'$\hat\eta$','$\tilde\eta$','$\hat\varphi$','$\tilde\varphi$'});
yticks([1 3 6 9 12]);
subplot(3,3,4)
boxplot([k1Hat,k1Tilde,k2Hat,k2Tilde,k3Hat,k3Tilde],'Symbol','');
set(gca,'TickLabelInterpreter','latex');
grid on
title('$T=60$','interpreter','latex');
ylim([0.95 12.3]);
xticklabels({'$\hat k_1$','$\tilde k_1$','$\hat k_2$','$\tilde k_2$',...
             '$\hat k_3$','$\tilde k_3$'});
yticks([1 3 6 9 12]);

%%% T = 120
T = 120;
a1 = (T-N-2)/N;
a2 = a1*(T-N-1)*(T-N-4)/(T-2);
a3 = a2/T;
M = eye(T)-(ones(T,1)*ones(1,T))./T;
rho = N/T;
NumberWindow = Ttotal-T;
EtaHat = zeros(NumberWindow,1);
VarphiHat = zeros(NumberWindow,1);
k1Hat = zeros(NumberWindow,1);
k2Hat = zeros(NumberWindow,1);
k3Hat = zeros(NumberWindow,1);
EtaTilde = zeros(NumberWindow,1);
VarphiTilde = zeros(NumberWindow,1);
k1Tilde = zeros(NumberWindow,1);
k2Tilde = zeros(NumberWindow,1);
k3Tilde = zeros(NumberWindow,1);
for j = 1:NumberWindow
    if mod(NumberWindow-j,100) == 0
        WindowsLeftT120 = NumberWindow-j
    end    
    X = R(j:j+T-1,:);
    %(\eta,\varphi) Student t
    [~,~,nuhat] = mlstudent(X);
    nuhat = min(nuhat,1000);
    y = @(x) (nuhat-2).*rho.*x./(2*(1-rho));
    fun = @(x) y(x).*exp(y(x)).*Expfun(nuhat/2,y(x))-rho;
    EtaHat(j,1) = fzero(fun,0.5*(1+nuhat/(nuhat-2)));
    VarphiHat(j,1) = 2*EtaHat(j,1)^2*(1-rho)/(nuhat-EtaHat(j,1)*(nuhat-2));
    %(\eta,\varphi) Elliptical
    muhat = mean(X);
    sigmahat = cov(X,1);
    tauhat = sum((X-muhat).^2,2);
    tauhat = tauhat./mean(tauhat);
    fun = @(x) sum(1./(T-N+x.*tauhat.*N))-1;
    EtaTilde(j,1) = fzero(fun,[1,10]);
    VarphiTilde(j,1) = (1-rho)./(1/EtaTilde(j,1)^2-sum((N.*tauhat.^2)./...
                       (T-N+N.*EtaTilde(j,1).*tauhat).^2));
    %(k1,k2,k3) Student t
    Expk1vec = zeros(nobs,1);
    Expk2vec = zeros(nobs,1);
    Expk3vec = zeros(nobs,1);
    for n = 1:nobs
         lambdahat = diag(sqrt((nuhat-2)./chi2rnd(nuhat,[T,1])));
         B = lambdahat*M*lambdahat;
         oneslambda = lambdahat*ones(T,1);
         Y = randn(T,N);
         Mat1 = inv(Y'*B*Y);
         Mat2 = Mat1^2;
         Expk1vec(n,1) = trace(Mat1);
         Expk2vec(n,1) = trace(Mat2);
         Expk3vec(n,1) = oneslambda'*Y*Mat2*Y'*oneslambda;
    end
    k1Hat(j,1) = a1*mean(Expk1vec);
    k2Hat(j,1) = a2*mean(Expk2vec);
    k3Hat(j,1) = a3*mean(Expk3vec);
    %(k1,k2,k3) Elliptical
    lambdahat = diag(sqrt(tauhat));
    B = lambdahat*M*lambdahat;
    oneslambda = lambdahat*ones(T,1);
    Expk1vec = zeros(nobs,1);
    Expk2vec = zeros(nobs,1);
    Expk3vec = zeros(nobs,1);
    for n = 1:nobs
         Y = randn(T,N);
         Mat1 = inv(Y'*B*Y);
         Mat2 = Mat1^2;
         Expk1vec(n,1) = trace(Mat1);
         Expk2vec(n,1) = trace(Mat2);
         Expk3vec(n,1) = oneslambda'*Y*Mat2*Y'*oneslambda;
    end
    k1Tilde(j,1) = a1*mean(Expk1vec);
    k2Tilde(j,1) = a2*mean(Expk2vec);
    k3Tilde(j,1) = a3*mean(Expk3vec);
end
subplot(3,3,2)
boxplot([EtaHat,EtaTilde,VarphiHat,VarphiTilde],'Symbol','');
set(gca,'TickLabelInterpreter','latex');
grid on
title('$T=120$','interpreter','latex');
ylim([0.95 5]);
xticklabels({'$\hat\eta$','$\tilde\eta$','$\hat\varphi$','$\tilde\varphi$'});
yticks([1 2 3 4 5]);
subplot(3,3,5)
boxplot([k1Hat,k1Tilde,k2Hat,k2Tilde,k3Hat,k3Tilde],'Symbol','');
set(gca,'TickLabelInterpreter','latex');
grid on
title('$T=120$','interpreter','latex');
ylim([0.95 5.2]);
xticklabels({'$\hat k_1$','$\tilde k_1$','$\hat k_2$','$\tilde k_2$',...
             '$\hat k_3$','$\tilde k_3$'});
yticks([1 2 3 4 5]);

%%% T = 240
T = 240;
a1 = (T-N-2)/N;
a2 = a1*(T-N-1)*(T-N-4)/(T-2);
a3 = a2/T;
M = eye(T)-(ones(T,1)*ones(1,T))./T;
rho = N/T;
NumberWindow = Ttotal-T;
EtaHat = zeros(NumberWindow,1);
VarphiHat = zeros(NumberWindow,1);
k1Hat = zeros(NumberWindow,1);
k2Hat = zeros(NumberWindow,1);
k3Hat = zeros(NumberWindow,1);
EtaTilde = zeros(NumberWindow,1);
VarphiTilde = zeros(NumberWindow,1);
k1Tilde = zeros(NumberWindow,1);
k2Tilde = zeros(NumberWindow,1);
k3Tilde = zeros(NumberWindow,1);
for j = 1:NumberWindow
    if mod(NumberWindow-j,100) == 0
        WindowsLeftT240 = NumberWindow-j
    end    
    X = R(j:j+T-1,:);
    %(\eta,\varphi) Student t
    [~,~,nuhat] = mlstudent(X);
    nuhat = min(nuhat,1000);
    y = @(x) (nuhat-2).*rho.*x./(2*(1-rho));
    fun = @(x) y(x).*exp(y(x)).*Expfun(nuhat/2,y(x))-rho;
    EtaHat(j,1) = fzero(fun,0.5*(1+nuhat/(nuhat-2)));
    VarphiHat(j,1) = 2*EtaHat(j,1)^2*(1-rho)/(nuhat-EtaHat(j,1)*(nuhat-2));
    %(\eta,\varphi) Elliptical
    muhat = mean(X);
    sigmahat = cov(X,1);
    tauhat = sum((X-muhat).^2,2);
    tauhat = tauhat./mean(tauhat);
    fun = @(x) sum(1./(T-N+x.*tauhat.*N))-1;
    EtaTilde(j,1) = fzero(fun,[1,10]);
    VarphiTilde(j,1) = (1-rho)./(1/EtaTilde(j,1)^2-sum((N.*tauhat.^2)./...
                       (T-N+N.*EtaTilde(j,1).*tauhat).^2));
    %(k1,k2,k3) Student t
    Expk1vec = zeros(nobs,1);
    Expk2vec = zeros(nobs,1);
    Expk3vec = zeros(nobs,1);
    for n = 1:nobs
         lambdahat = diag(sqrt((nuhat-2)./chi2rnd(nuhat,[T,1])));
         B = lambdahat*M*lambdahat;
         oneslambda = lambdahat*ones(T,1);
         Y = randn(T,N);
         Mat1 = inv(Y'*B*Y);
         Mat2 = Mat1^2;
         Expk1vec(n,1) = trace(Mat1);
         Expk2vec(n,1) = trace(Mat2);
         Expk3vec(n,1) = oneslambda'*Y*Mat2*Y'*oneslambda;
    end
    k1Hat(j,1) = a1*mean(Expk1vec);
    k2Hat(j,1) = a2*mean(Expk2vec);
    k3Hat(j,1) = a3*mean(Expk3vec);
    %(k1,k2,k3) Elliptical
    lambdahat = diag(sqrt(tauhat));
    B = lambdahat*M*lambdahat;
    oneslambda = lambdahat*ones(T,1);
    Expk1vec = zeros(nobs,1);
    Expk2vec = zeros(nobs,1);
    Expk3vec = zeros(nobs,1);
    for n = 1:nobs
         Y = randn(T,N);
         Mat1 = inv(Y'*B*Y);
         Mat2 = Mat1^2;
         Expk1vec(n,1) = trace(Mat1);
         Expk2vec(n,1) = trace(Mat2);
         Expk3vec(n,1) = oneslambda'*Y*Mat2*Y'*oneslambda;
    end
    k1Tilde(j,1) = a1*mean(Expk1vec);
    k2Tilde(j,1) = a2*mean(Expk2vec);
    k3Tilde(j,1) = a3*mean(Expk3vec);
end
subplot(3,3,3)
boxplot([EtaHat,EtaTilde,VarphiHat,VarphiTilde],'Symbol','');
set(gca,'TickLabelInterpreter','latex');
grid on
title('$T=240$','interpreter','latex');
ylim([0.95 2.5]);
xticklabels({'$\hat\eta$','$\tilde\eta$','$\hat\varphi$','$\tilde\varphi$'});
yticks([1 1.5 2 2.5]);
subplot(3,3,6)
boxplot([k1Hat,k1Tilde,k2Hat,k2Tilde,k3Hat,k3Tilde],'Symbol','');
set(gca,'TickLabelInterpreter','latex');
grid on
title('$T=240$','interpreter','latex');
ylim([0.95 2.5]);
xticklabels({'$\hat k_1$','$\tilde k_1$','$\hat k_2$','$\tilde k_2$',...
             '$\hat k_3$','$\tilde k_3$'});
yticks([1 1.5 2 2.5]);
fontsize(fig,16,"points");